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1.  INTRODUCTION 


Predicting  the  propagation  of  optical  frequency  beams  through  the  atmosphere 
depends  greatly  on  knowledge  of  the  detailed  state  of  the  atmosphere  along  the 
propagation  path.  This  statement  is  appropriate  regardless  of  the  altitude  at  which  the 
propagation  is  examined.  Two  atmospheric  characteristics  that  directly  impede  beam 
transmission  are  optical  scale  turbulence  and  the  presence  of  clouds  due  to  gravity  waves 
and  other  generation  processes.  It  has  also  been  observed  that  gravity  wave  instability 
processes  can  lead  to  local  turbulence  production.  However,  this  is  one  of  the  “-least 
quantified  aspects  of  gravity  wave  forcing  of  the  middle  atmosphere  at  present,”  as  noted 
by  Fritts  and  Alexander  in  a  recent  review  article  [1]. 

In  the  case  of  optical  turbulence,  considerable  experimental  effort  has  been  expended 
toward  expanding  the  data  base  of  temperature  fluctuations  in  the  atmosphere,  with 
accepted  average  profiles  having  been  obtained  for  selected  important  locales  as  noted  in 
Reference  2  which  provides  an  excellent  assessment  of  past  work  and  the  issues. 
Unfortunately,  measuring  directly  the  detailed  state  of  the  atmosphere  is  not  usually 
possible  or  practical  for  a  significant  range  of  conditions,  locations  and  times,  especially 
given  terrain  and  local  condition  influences.  This  situation  points  to  the  need  for 
simulation/modeling  tools  that  could  be  used  to  predict,  based  on  inputs  of  terrain  and 
local  conditions,  local  levels  of  optical  turbulence.  However,  the  presently  available 
atmospheric  modeling  tools  do  not  provide  a  means  of  simulating  directly  the  scales 
needed  to  assess  local  optical  turbulence,  due  to  resolution  and  other  issues. 

This  research  project  was  initiated  to  continue  development  of  a  high  resolution 
meso-scale  atmospheric  model,  based  on  the  well  developed  model  MM5  [3],  to  provide 
dynamic  resolution  of  selected  features  and  characteristics  of  the  atmosphere  simulation. 
The  goal  is  to  provide  a  very  high  resolution  prediction  of  optical  turbulence  parameters. 
The  remainder  of  this  report  will  review  the  code  development  and  improvements 
performed  during  the  first  year  of  the  contract  period  followed  by  a  review  of  the  results 
obtained  to  date.  Relevant  details  of  the  following  activity  will  be  included: 

1)  Installation  of  improved  boundary  condition  implementations  for  use  when  the 
adaptive  algorithm  is  used  for  the  entire  mesoscale  domain.  This  permits 
adaptation  to  be  used  instead  of  nesting  to  increase  resolution.  This 
implementation  will  also  provide  improved  one  way  nest  boundary  conditions. 

2)  Further  checkout  of  the  transformation  coding  was  performed. 

3)  A  new  method  was  installed  for  redistribution  of  the  solution  after  adaptation  of 
the  mesh. 

4)  Basic  development  of  the  atmospheric  version  of  the  NCSU  k-£  turbulence  model 
was  completed  during  this  period.  The  fully  developed  version  was  installed  in 
the  updated  code. 
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A  2-D  case  consisting  of  sounding  data  from  a  Colorado  windstorm  event  applied  to  a 
domain  containing  a  generic  mountain  [4]  was  used  for  code  development.  Many  other 
models  have  been  executed  for  this  case  and  the  strong  shear  coupled  with  the  generic 
mountain  topography  provides  conditions  usually  leading  to  gravity  waves  and 
turbulence.  Standard  confirmation  steps  for  the  adaptive  mesh  algorithm  were  performed 
using  this  case.  Currently,  3-D  runs  are  in  progress  to  compare  the  results  of  the  high- 
resolution  mesoscale  model  and  turbulence  model  with  data  and  other  atmospheric 
models. 

2.  DYNAMICALLY  RESOLVED  MESOSCALE  MODEL- 
MODIFICATIONS  TO  MM5 

The  code  structure  and  physical  modeling  in  MM5  [3]  is  retained.  A  3-D  grid  domain 
(nest)  is  imbedded  in  the  inner  MM5  nest  with  increased  mesh  node  density  and  with  the 
NCSU  dynamic  solution  adaptive  mesh  algorithm  DSAGA  [5-7]  applied  to  further 
increase  resolution  where  needed.  The  flow  in  this  imbedded  region  is  solved,  to  LES 
scales,  with  the  standard  non-hydrostatic  equations  transformed  to  a  general  time-varying 
structured  grid.  DSAGA  uses  r-refinement  adaptation  to  resolve  automatically  selected 
scales  and  features  in  the  solution.  The  goal  of  the  adaptation  is  to  provide  LES-scale 
resolution  for  the  developing  turbulence.  This  capability  will  also  allow  resolution  of  the 
dynamic  processes  from  gravity  wave  generation  to  break-down.  This  adaptive  algorithm 
is  developed  sufficiently  such  that  any  criteria  or  linear  combination  of  criteria  may  be 
used  to  promote  grid  clustering,  thereby  insuring  resolution  of  the  features  from  which 
the  criteria  are  derived.  Within  the  dynamically  resolved  nest,  an  extension  of  k-  C, 
(enstrophy,  or  variance  of  vorticity)  turbulence  model  [8,  9]  is  used  to  provide  LES  sub¬ 
grid  scale  turbulence  modeling.  This  model  includes  fundamental  flow  physics  and  has 
proved  to  be  superior  to  standard  models  when  used  in  modeling  turbulence. 

The  outer  MM5  and  imbedded  high-resolution  fields  are  executed  alternately.  After 
each  iteration  on  the  outer  MM5  domain,  the  boundary  values  and  tendencies  for  the 
imbedded  mesh  are  interpolated  from  the  outer  MM5  domain.  Temporal  accuracy  is 
preserved  by  advancing  the  solution  in  the  imbedded  domain  to  the  same  time  level  as  the 
outer  domain.  This  procedure  is  the  same  as  MM5,  but  the  time  step  in  the  imbedded 
domain  may  vary  due  to  the  grid  adaptation. 

3.  TRANSFORMATION  OF  GOVERNING  EQUATIONS 

In  order  to  install  the  dynamic  adaptive  mesh  algorithm  in  MM5,  the  non-hydrostatic 
governing  equations  were  transformed  to  a  time-varying  curvilinear  coordinate  system 
that  is  coincident  with  either  of  the  outer  domain,  a  nest  boundary,  or  a  portion  of  the 
standard  nest.  The  transformation  is  Galilean  invariant,  implying  mathematically  that  the 
solution  to  the  governing  equations  is  not  changed  by  the  mesh  adjustment  to  improve 
solution  resolution. 
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The  equations,  initially  defined  in  the  x,  y,  cr  coordinate  system,  are  transformed  using 
the  chain-rule  in  all  three  dimensions  to  a  uniform  computational  coordinate  system, 
according  to: 

T  =  t 

4,=^(x,y,z,a) 

where  cr  is  the  non-dimensional  pressure  coordinate.  The  resulting  equations  read 

dU  mat  2  dE  34,  ,  n  2  5E  34,  ,  5G  34,  _c 

- h - b  TYl - h  171 - 1 - —  o 

dt  d4j  dt  d^j  dx  5£,  dy  dd,,  do 

where 

E  =  [p*  p',p*u,p*v,p*w,p*Tf 

E  =  —  [p*p',p*u,p*v,p*w,p*T]r 
m 

F  =  —  \p*  p',p*u,p*  v,p*w,p*Tj 
m 

G  =  &[p*  p' ,p*  u,p*  v,p*  w,p*Tj 


m  is  the  map  scale,  p*  the  reference  pressure,  and  p  the  pressure  perturbation,  T  the 
temperature,  and  u,  v  and  w  are  velocity  components  in  the  x-,  y-  and  z-direction, 
respectively.  All  other  terms,  such  as  pressure  gradient,  Coriolis  force,  and  gravity  terms, 
are  included  in  S,  c.f.  Reference  3  for  more  detail.  The  above  equations  are  discretized  in 
the  Arakawa-B  [10]  type  staggered  grid,  using  the  same  finite  difference  stencils  as 
MM5,  e.g.,  the  stencils  used  in  the  x-direction  in  MM5  are  applied  to  direction  here. 
These  equations  are  also  solved  using  the  leap-frog  scheme.  In  order  to  obtain  accurate 
discretization  in  the  curvilinear  staggered  grid,  three  sets  of  metric  derivatives  are 
calculated  to  be  consistent  to  the  differencing  of  flow  variables.  The  variables  and  metric 
derivatives  are  defined  at  three  different  locations  as  shown  in  Figure  1,  i.e.,  the  cell 
center  (cross,  for  p'  and  T),  the  center  of  cell  edges  (dot,  for  u  and  v)  and  the  center  of 
43  =  const  cell  surface  ( A ,  for  w  and  cr ). 


4.  MODIFICATION  OF  THE  SEMI-IMPLICIT  SCHEME 


In  MM5,  in  order  to  remove  the  limitation  on  the  time  step  due  to  small  mesh  spacing 
in  the  vertical  direction,  the  following  two  coupled  equations  for  w  and  p'  are  solved 
implicitly: 

aw  p0g  dp1 !  g  p'  _s 
dt  pp*  da  y  P 


dt 


Pogyp  dw 

p*  d(j 


~PQgw  = 


S 


p' 
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This  results  in  a  tridiagonal  system  along  the  a  -direction  for  W  in  the  uniform  mesh, 
which  can  be  solved  directly.  But  in  the  adaptive  mesh,  the  cr  variation  is  transformed, 
therefore,  this  set  must  be  solved  iteratively.  The  iteration  scheme  chosen  is  as  follows: 

1)  Initialize  as: 


2)  Then  solve  iteratively: 

w('+l)  —  w1  Ppgffo'Y  } 
A/  PP*[d^j 


%3,e  + 


g  P 


■(<+)) 


y  p 


— s„  + 


Pag 


PP 


'w'v' 

&,j 


+ 


p'{M)-p'‘  p0gyp 


At 


P* 


^-Pa^=Sp^ 


Pagyp 


\PtxJ 


\d$iJ 

£l,cr  + 


\dZ2  J 


3)  When  converged. 


p'M  =  p*+l\w‘+l=wm) 


5.  ADAPTIVE  GRID  ALGORITHM 

In  the  current  implementation  of  the  adaptive  algorithm,  the  integration  proceeds  in  two 
steps;  the  first  is  the  integration  of  the  governing  equations  on  the  mesh  adapted  to  the 
solution  at  the  previous  time  step.  After  this  integration  is  completed,  the  mesh  is  adapted 
to  this  new  solution  and  the  invariant  solution  is  redistributed  over  the  new  mesh.  The 
weight  function  that  controls  the  adaptation  is  based  on  the  magnitude  of  local  vorticity 
and  is  processed  and  smoothed  to  influence  the  relative  adaptation  of  the  grid  and  the 
smoothness  of  the  grid  metric  derivative  distribution.  The  redistribution  step  is  equivalent 
to  an  advection  of  the  invariant  solution  and  is  thereby  important  to  the  accuracy  of  the 
integration.  In  order  to  improve  the  accuracy  of  this  step,  a  fifth  order  weighted 
essentially  non-oscillatory  (WENO)  scheme  [11]  was  installed.  As  confirmation  of  this 
advection  algorithm  and  as  final  confirmation  of  the  metric  derivatives  of  the 
transformation,  solutions  were  obtained  using  the  unmodified  MM5  and  using  the 
transformed  equations  and  adaptive  algorithm.  Overlay  of  the  solution  contours  should 
provide  confirmation  of  the  Galilean  invariance  of  the  adaptive  solution  and  of  the 
metrics.  As  shown  in  Fig.  2,  the  adaptive  procedure  recovers  the  same  solution  as  that  of 
the  fixed  mesh  except  where  the  increased  resolution  of  the  adaptive  mesh  improves 
solution  detail. 

6.  PROGRESS  ON  TURBULENCE  MODELING 

During  this  reporting  period,  the  turbulence  modeling  was  completed  and  runs  were 
started.  The  k- £  model  has  been  expanded  to  include  equations  for  the  enthalpy  variance 
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Figure  1 .  Arakawa  B  grid  showing  variable  location. 


Figure  2.  Comparison  of  potential  temperature  results  (Colorado  windstorm  conditions 
[2]);  dotted  line:  solution  on  the  uniform  mesh;  dashed  line:  solution  on  the 
adaptive  mesh  after  applying  the  WENO  redistribution. 
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and  its  dissipation  rate.  The  modeled  equations  are  included  in  Appendix  A  in  complete 
form. 


A  hybrid  LES/RANS  version  of  k model  is  given  by 
D(pk ) 


Dt 


=  V« 


(  \ 

E+£l 

Vifc 

K  3  ak  J 

Mi  g  S0 


+  r  ::S 

Pr,  6  dz 


(i  -  r)(-^  vp .  v/>  +  c,  ^ + /<)  -  rc,p 

ck  P  Tp  A 

The  blending  function  (T),  that  blends  the  LES  and  RANS  and  controls  the  regions 
where  LES  or  RANS  is  dominant.  In  Ref.  9,  the  following  blending  was  developed  for 
wall-bounded  flows: 

f  i  Y 

Ta  =  tanh  — — 


where  A  is  the  minimum  mesh  spacing,  ls  is  the  turbulence  length  scale,  ai  the  model 
constant  used  to  scale  the  influence  of  the  minimum  mesh  size.  Although  this  constant 
was  originally  modeled  using  aerodynamic  flows,  it  did  not  function  properly  when  used 
for  atmospheric  turbulence.  The  LES  subgrid  model  has  been  used  in  the  test  cases  in  this 
report. 


7.  RESULTS  AND  DISCUSSION 


7.1.  2-D  Case 


A  two-dimensional  case,  based  on  an  actual  frontal  event  and  subsequent  windstorm 
downslope  of  the  Colorado  front  range,  is  used  to  investigate  the  ability  of  DSAGA  to 
resolve  shear  layers,  gravity  waves  and  their  breakdown.  This  case  is  the  same  as  that 
used  in  Reference  4  for  evaluating  the  wave  breaking  prediction  for  various  models  by 
calculating  the  response  to  a  generic  mountain  in  the  presence  of  an  actual  upwind 
sounding.  For  this  2-D  case,  the  adaptive  flow  solver  is  executed  on  the  full  domain.  The 
mountain  profile  is  prescribed  by  a  witch  of  Agnesi  curve 


zs(x)  = 


ha 2 


x2  +  a 
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where  h  =  2  km  is  the  mountain  height  and  a  —  10  km  is  the  mountain  half  width,  as 
shown  in  Figure  3.  In  order  to  compare  with  those  results  in  Reference  4,  the  initial 
conditions  are  the  same  as  theirs  -  horizontally  homogeneous  with  the  initial  flow  profile 
provided  from  the  upstream  sounding  data  of  11  January  1972,  Grand  Junction, 
Colorado,  as  shown  in  Figure  4.  The  same  uniform  grid  is  used  to  initialize  the 
computation,  which  has  221  x  126  nodes  with  Az  ~  200  m  and  Ax  ~  1  km.  Free  slip 
boundary  conditions  are  applied  to  the  lower  boundary.  Orlanski  radiation  boundary 
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conditions  [12]  are  applied  to  both  inflow  and  outflow  boundaries.  A  5  km  sponge  layer 
is  set  at  the  top  of  the  domain  as  a  boundary  condition  to  prevent  the  reflection  of  waves 
from  the  top  boundary.  A  fourth  order  damping  term  is  added  to  momentum  and 
temperature  equations  to  stabilize  the  flow  solver. 

In  the  following  results,  the  grid  adaptation  is  started  from  t  =  0.  The  grid  is 
preadapted  from  the  uniform  mesh  by  100  adaptation  cycles  in  order  to  resolve  the  initial 
flow  field.  The  4  grid  adaptation  cycles  were  performed  after  every  4  flow-solver 
iterations.  One  grid  adaptation  cycle  is  defined  as  the  sequence  of  weight-function 
calculation,  grid  point  relocation  and  solution  redistribution.  Two  sets  of  results  are 
obtained:  one  without  turbulence  modeling,  the  other  with  our  hybrid  RANS/LES  model. 
For  the  former,  the  inherent  numerical  damping  in  the  adaptive  solver  provides 
dissipation  mechanism  for  turbulence,  which  can  be  considered  as  a  Monotonically 
Integrated  LES  (MILES)  simulation  [13].  For  the  hybrid  RANS/LES  simulation,  some 
tests  showed  that  the  hybrid  scheme  cannot  switch  from  RANS  to  LES  and  turbulent 
kinetic  energy  (TKE).  Figure  5  shows  the  potential  temperature  contours  (0),  or  the 
isotropes,  for  the  MILES  simulation  at  /  =  5000  sec.  Compared  to  the  uniform  mesh 
solution,  the  isotropes  on  the  adaptive  mesh  clearly  show  the  dynamic  process  of  gravity 
wave  overturning  and  break-down  to  turbulent  eddies  at  1 8  <  z  <  20  km  and  13  <  z  <  16 
km.  Note  that  the  breakdown  is  occurring  in  two  distinct  shear  layers.  The  coherent 
turbulent  structure  at  upper  level  atmosphere  is  essentially  absent  from  the  uniform  mesh 
solutions  because  the  1  km  mesh  spacing  is  not  able  to  resolve  the  1 .5  km  long  eddies  at  z 
~  20  km.  As  shown  in  Figure  6,  the  wave  break-down  is  closely  correlated  with  the 
regions  of  high  vorticity  (thus  the  choice  of  vorticity  for  the  weight  function).  As  a  result 
of  this,  the  grid  nodes  cluster  around  the  wave  breaking  region  and  provide  local 
resolution  of  Ax  ~  300  m  and  Az  ~  50  m  for  resolving  the  small  eddies.  The  potential 
temperature  contours  at  t  =  3  h  for  RANS/LES  on  the  adaptive  mesh  are  shown  in  Figure 
7.  The  adaptive  mesh  at  that  moment  is  compared  with  the  uniform  mesh  in  Figure  8. 
Two  bands  of  highly  refined  mesh  at  z  ~  13  km  and  20  km  above  the  mountain  provide 
better  resolution  on  small  eddies  within  those  regions.  Another  region  of  refined  grid  is 
aligned  with  the  vortex  street  above  the  lee  side  of  the  mountain  gives  better  resolution  of 
lower-tropospheric  eddies  in  front  of  the  hydraulic  jump.  The  third  high  resolution  region 
is  located  on  the  top  of  the  hydraulic  jump  due  to  a  large  size  eddy  located  there.  Figure  9 
shows  the  contour  of  aircraft  lift  change  response  (to  vertical  gusts/velocity  components), 
which  can  be  used  as  an  index  to  assess  the  turbulence  level.  The  lift  change  is  defined 
as: 

AZ,  =  V  c 
L  230 (m/s)  La 

where  the  AL/Lis  the  relative  change  of  lift,  V  the  vertical  velocity  component  and 
CLa  =  6  is  a  lift  coefficient.  As  can  be  seen  from  the  AZ,  /  L  contour,  the  lift  change  can 

be  as  large  as  1  g,  which  has  an  adverse  effect  on  aviation  operation  and  safety  in  this 
region. 

Figure  10  gives  a  profile  of  C2n  at  x  =  30  km,  t  =  3  h,  which  is  calculated  using 
Dewan’s  formula  [14].  When  compared  to  Figure  8,  there  is  some  correlation  between 
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Figure  5.  Comparison  of  potential  temperature  contours  on  fixed  and  adapted  mesh, 
LES-RANS. 


weight  function  adaptive  mesh 
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Figure  6.  Typical  response  of  adaptive  mesh  to  weight  function,  LES-RANS. 
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Figure  7.  Potential  temperature  contours  at  t  =  3  h. 


Figure  8.  Detail  of  adapted  computational  mesh  at  the  t=3  h  solution  time. 
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Figure  9.  Lift  change  (A L/  L)  contours 


Figure  10.  Cj;  profile  at  x  =  30  km,  t  =  3  h. 
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high  Cl  regions  and  refined  mesh  regions.  Therefore,  these  high  C„2  regions  are  due 
mainly  to  high  shear  and  wave-breaking. 


7.2.  3-D  Case 


In  this  case,  the  adaptive  algorithm  is  integrated  into  MM5  code  to  simulate  the  SOC 
(Storm  of  the  Century)  case,  a  tutorial  case  for  MM5.  Figure  1 1  presents  the  layout  of  the 
two-level  nested  grid  for  this  simulation  with  MM5  36employed  in  both  levels.  Its 
solution  in  level  2  provides  the  boundary  conditions  for  the  current  adaptive  flow  solver. 
The  initial  grid  for  the  adaptive  solver  is  the  same  as  the  level  2  grid  with  49x52x23 
nodes.  For  this  case,  all  Planetary  Boundary  Layer  options  in  MM5  are  turned  off  and  the 
free  slip  boundary  condition  is  used.  The  MILES  scheme  is  used  in  the  imbedded  region. 
To  advance  the  solution  to  the  same  level  as  the  outer  MM5  solution,  the  time  step  in  the 
imbedded  adaptive  mesh  is  determined  by: 


At  =  At 


MM  5 


/  nint( 


At 


MM  5 


min(A thc) 


0.5) 


where  A tMMi  is  the  time  step  on  the  outer  MM5  domain,  A the  is  the  local  time  step,  nint 
is  a  Fortran  function  for  nearest  integer.  Figure  12  shows  the  evolution  of  the  grid  in  a 
£,=  const  plane.  The  contours  for  the  magnitude  of  vorticity  are  also  shown  in  the  same 
figure.  As  can  be  seen  from  the  grid  node  distribution,  the  high  shear  layers  near  the 
surface  and  in  the  upper  atmosphere  are  dynamically  captured,  i.e.,  the  high  resolution 
grid  always  adapts  to  the  region  with  large  vorticity  (shear)  and  moves  with  it.  Figure  13 
shows  the  grid  and  contours  in  £3  =  const  plane  next  to  the  surface.  Again,  the  high 

vorticity  region  is  well  resolved.  The  minimum  mesh  spacing  in  the  region  shown  in 
Figure  13  is  about  6  km,  which  is  about  one  fifth  of  the  initial  mesh  spacing  (30  km). 
These  two  figures  clearly  demonstrate  that  the  current  algorithm  can  resolve  high  velocity 
gradients  in  both  vertical  and  horizontal  directions. 


8.  CONCLUSIONS 


Significant  progress  was  made  during  this  reporting  period  towards  the  ultimate  goal 
of  a  fully  functioning  mesoscale  model  for  predicting  the  likelihood  of  optical  turbulence 
under  specified  conditions  and  topography.  Incorporation  of  the  k-C,  turbulence  model 
into  the  entire  modeling  system  was  completed  and  checkout  of  the  model  blending 
algorithm  indicates  some  differences  in  scale  as  might  be  expected.  Reproduction  of  the 
fixed  mesh  results  by  the  adaptive  mesh  algorithm  was  confirmed.  The  results  show  that 
the  vorticity  is  a  good  indicator  for  gravity  wave-breaking.  The  entire  wave-breaking 
process  in  the  2-D  case  is  well-captured  and  resolved  in  adaptive  mesh. 
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Figure  1 1 .  The  layout  of  two-level  nested  grid  and  terrain  contours. 
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Figure  12.  Grid  and  the  vorticity  magnitude  contours  in  ^  =  const,  surface. 
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Figure  13.  The  vorticity  magnitude  contours  and  horizontal  velocity  vectors  in  =  const 
plane  at  t  =  720  min. 
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Appendix. 

Turbulence  Model  Equations 
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Enthalpy  Variance  Equation: 
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where 
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Dissipation  Rate  Equation  for  Enthalpy  Variance: 
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